#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "choletsky.h"
#include <time.h>
#include <cublas.h>

int N = 1000;


void printFirstPart(int n, int r, double ea, double eb)
{
	printf("\n     P R O B L E M:\n");
	printf("         solving of the  linear algebraic system\n");
	printf("         with a symmetic positive defined matrix \n");
	printf("\n     D a t a :\n");
	printf("         - number of matrix's rows                  = %d",n);
	printf("\n");
	printf("         - number of matrix's columns               = %d",n);
	printf("\n");
	printf("         - number of the right-hand side             ");
	printf("\n");
	printf("                   of the systems                   = %d",r);
	printf("\n");
	printf("         - maximum relative error              ");
	printf("\n");
	printf("                   of the matrix elements           = %.5e",ea);
	printf("\n");
	printf("         - maximum relative error              ");
	printf("\n");
	printf("           of elements of the right-hand sides      = %.5e",eb);
	printf("\n");
	printf("\n");
	printf("    P r o c e s s   o f   i n v e s t i g a t i n g   a n d    s o l v i n g"); 
	printf("\n");
	printf("\n     M E T H O D:   \n");
	printf("         Choletsky decomposition   \n \n");
	printf(" R E S U L T S :\n");

}

void save_result(char * file_name, double *res)
{
	FILE *res_file;
	if ( (res_file =fopen(file_name,"w"))==NULL)
	{
		printf("can not create result file");
	}
	else
	{
		
		for(int i =0; i <N;i++)
		{
			fprintf(res_file,"%f\n",res[i]);	
		}
		fclose(res_file);
	}
}

int main(int argc, char** argv)
{   
    double* A = NULL;
    double* B = NULL;
    double* X = NULL;
    int* IPIV = NULL;
	double *E;
	double *E1;
	double *cond;
	double ea, eb;
    int i, j;
	clock_t start_time, end_time;
   
    if(argc >=  2)
	{
		printf(argv[1]);
		printf("\n!!");
		N = atoi(argv[1]);
	}
    A = (double*)malloc(N*N*sizeof(double));
    B = (double*)malloc(N*sizeof(double));
    X = (double*)malloc(N*sizeof(double));
	E = (double *)malloc(sizeof(double));
	E1 = (double *)malloc(sizeof(double));
	cond = (double *)malloc(sizeof(double));

	ea = 0;
	eb = 0;
    
    for(i = 0; i < N; ++i)
	{    
		for(j = 0; j < N; ++j)
        {       
            A[i*N+j] = i>j? (N- i):(N - j);
        }
	}

    for(i = 0; i < N; ++i)
    {
        B[i] = 0;
    }
    B[1] = 1.0f;

    memcpy(X, B, N*sizeof(double));
   
	printFirstPart( N, 1, ea, eb);

	int ret_value;	
	start_time = clock();

	ret_value = LPPSADB_cu(N, A, X, E, E1, cond, ea, eb);

	end_time = clock();


	if (ret_value == 0)
	{
		save_result("result.out", X);
		printf(" \n SOLUTION IS OBTAINED IN FILE result.out \n");

		int q1=5,q2=N-5;

		printf(" \n \n Solution (first and last 5 components) for right part number\n \n");
		for(i=0;i<q1;i++)
			printf("%25.16e ",X[i]);
		printf(" \n \n");
		for(i=q2;i<N;i++)
			printf("%25.16e ",X[i]);
		printf("\n");

		printf("\n       E s t i m a t e s : ");
		printf("\n       - inherited error in the solution : %.7e",*E );
		printf("\n       - computational error in the solution: %.7e",*E1 );
		printf("\n       - matrix condition number: %.7e",*cond );
		printf("\nTime : %.5f\n",(float)(end_time - start_time)/CLOCKS_PER_SEC);

	}
	else if(ret_value == -1)
	{
		printf("        !!! THE MATRIX IS MACHINE-SINGULAR\n");
		printf("\nTime : %.5f\n",(float)(end_time - start_time)/CLOCKS_PER_SEC);
	}
	else if(ret_value == -2)
	{
		printf("\n  !!!  W A R N I N G : \n");
		printf("\n");
		printf("         the solution's reliability isn't guaranteed\n");
		printf("         since the matrix may turn out to be singular \n");
		printf("         within the accracy of its elements' specification.\n");
		printf("         Esimates of the solution are not evaluated.\n");

		int q1=5,q2=N-5;

		printf(" \n \n Solution (first and last 5 components) for right part number\n \n");
		for(i=0;i<q1;i++)
			printf("%25.16e ",X[i]);
		printf(" \n \n");
		for(i=q2;i<N;i++)
			printf("%25.16e ",X[i]);
		printf("\n");
		printf("\nTime : %.5f\n",(float)(end_time - start_time)/CLOCKS_PER_SEC);
	}
	else if(ret_value == -4)
	{
		printf("        !!! THE MATRIX IS NOT POSITIVE DEFINED\n");
		printf("\nTime : %.5f\n",(float)(end_time - start_time)/CLOCKS_PER_SEC);
	}

    free(A);
    free(B);
	free(X);
	free(E);
	free(E1);
	free(cond);

}

